function [gi git] = h2ginverse(params, g, gt)

  gi  = zeros(params.Nx, params.Ny, params.Nz, 4, 4);
  git = zeros(params.Nx, params.Ny, params.Nz, 4, 4);

  detg = g(:,:,:,1,3).^2.*g(:,:,:,2,4).^2+(-2).*g(:,:,:,1,4).*g(:,:,:,1,3).*g(:,:,:,2,4).*g(:,:,:,2,3)+g(:,:,:,1,4).^2.*g(:,:,:,2,3).^2+(-2).*g(:,:,:,1,2).*g(:,:,:,1,3).*g(:,:,:,2,4).*g(:,:,:,3,4)+2.*g(:,:,:,1,4).*g(:,:,:,1,3).*g(:,:,:,2,2).*g(:,:,:,3,4)+(-2).*g(:,:,:,1,4).*g(:,:,:,1,2).*g(:,:,:,2,3).*g(:,:,:,3,4)+2.*g(:,:,:,1,1).*g(:,:,:,2,4).*g(:,:,:,2,3).*g(:,:,:,3,4)+g(:,:,:,1,2).^2.*g(:,:,:,3,4).^2+(-1).*g(:,:,:,1,1).*g(:,:,:,2,2).*g(:,:,:,3,4).^2+2.*g(:,:,:,1,4).*g(:,:,:,1,2).*g(:,:,:,2,4).*g(:,:,:,3,3)+(-1).*g(:,:,:,1,1).*g(:,:,:,2,4).^2.*g(:,:,:,3,3)+(-1).*g(:,:,:,1,4).^2.*g(:,:,:,2,2).*g(:,:,:,3,3);

  gi(:,:,:,1,1) = 2.*g(:,:,:,2,4).*g(:,:,:,2,3).*g(:,:,:,3,4)+(-1).*g(:,:,:,2,2).*g(:,:,:,3,4).^2+(-1).*g(:,:,:,2,4).^2.*g(:,:,:,3,3);
  gi(:,:,:,1,2) = (-1).*g(:,:,:,1,3).*g(:,:,:,2,4).*g(:,:,:,3,4)+(-1).*g(:,:,:,1,4).*g(:,:,:,2,3).*g(:,:,:,3,4)+g(:,:,:,1,2).*g(:,:,:,3,4).^2+g(:,:,:,1,4).*g(:,:,:,2,4).*g(:,:,:,3,3);
  gi(:,:,:,1,3) = g(:,:,:,1,3).*g(:,:,:,2,4).^2+(-1).*g(:,:,:,1,4).*g(:,:,:,2,4).*g(:,:,:,2,3)+(-1).*g(:,:,:,1,2).*g(:,:,:,2,4).*g(:,:,:,3,4)+g(:,:,:,1,4).*g(:,:,:,2,2).*g(:,:,:,3,4);
  gi(:,:,:,1,4) = (-1).*g(:,:,:,1,3).*g(:,:,:,2,4).*g(:,:,:,2,3)+g(:,:,:,1,4).*g(:,:,:,2,3).^2+g(:,:,:,1,3).*g(:,:,:,2,2).*g(:,:,:,3,4)+(-1).*g(:,:,:,1,2).*g(:,:,:,2,3).*g(:,:,:,3,4)+g(:,:,:,1,2).*g(:,:,:,2,4).*g(:,:,:,3,3)+(-1).*g(:,:,:,1,4).*g(:,:,:,2,2).*g(:,:,:,3,3);
  gi(:,:,:,2,2) = 2.*g(:,:,:,1,4).*g(:,:,:,1,3).*g(:,:,:,3,4)+(-1).*g(:,:,:,1,1).*g(:,:,:,3,4).^2+(-1).*g(:,:,:,1,4).^2.*g(:,:,:,3,3);
  gi(:,:,:,2,3) = (-1).*g(:,:,:,1,4).*g(:,:,:,1,3).*g(:,:,:,2,4)+g(:,:,:,1,4).^2.*g(:,:,:,2,3)+(-1).*g(:,:,:,1,4).*g(:,:,:,1,2).*g(:,:,:,3,4)+g(:,:,:,1,1).*g(:,:,:,2,4).*g(:,:,:,3,4);
  gi(:,:,:,2,4) = g(:,:,:,1,3).^2.*g(:,:,:,2,4)+(-1).*g(:,:,:,1,4).*g(:,:,:,1,3).*g(:,:,:,2,3)+(-1).*g(:,:,:,1,2).*g(:,:,:,1,3).*g(:,:,:,3,4)+g(:,:,:,1,1).*g(:,:,:,2,3).*g(:,:,:,3,4)+g(:,:,:,1,4).*g(:,:,:,1,2).*g(:,:,:,3,3)+(-1).*g(:,:,:,1,1).*g(:,:,:,2,4).*g(:,:,:,3,3);
  gi(:,:,:,3,3) = 2.*g(:,:,:,1,4).*g(:,:,:,1,2).*g(:,:,:,2,4)+(-1).*g(:,:,:,1,1).*g(:,:,:,2,4).^2+(-1).*g(:,:,:,1,4).^2.*g(:,:,:,2,2);
  gi(:,:,:,3,4) = (-1).*g(:,:,:,1,2).*g(:,:,:,1,3).*g(:,:,:,2,4)+g(:,:,:,1,4).*g(:,:,:,1,3).*g(:,:,:,2,2)+(-1).*g(:,:,:,1,4).*g(:,:,:,1,2).*g(:,:,:,2,3)+g(:,:,:,1,1).*g(:,:,:,2,4).*g(:,:,:,2,3)+g(:,:,:,1,2).^2.*g(:,:,:,3,4)+(-1).*g(:,:,:,1,1).*g(:,:,:,2,2).*g(:,:,:,3,4);
  gi(:,:,:,4,4) = (-1).*g(:,:,:,1,3).^2.*g(:,:,:,2,2)+2.*g(:,:,:,1,2).*g(:,:,:,1,3).*g(:,:,:,2,3)+(-1).*g(:,:,:,1,1).*g(:,:,:,2,3).^2+(-1).*g(:,:,:,1,2).^2.*g(:,:,:,3,3)+g(:,:,:,1,1).*g(:,:,:,2,2).*g(:,:,:,3,3);



  % notice still needs to divide gi by det AND symmetrize





  % divide by the determinant
  for i=1:4
    for j=i:4
      gi(:,:,:,i,j) = gi(:,:,:,i,j)./detg;
    end
  end




%  keyboard()
  % fill the symm part
  for i=1:4
    for j=i+1:4
      gi(:,:,:,j,i) = gi(:,:,:,i,j);
    end
  end


  tmp = zeros(params.Nx, params.Ny, params.Nz, 4, 4);

  % git = gi gt gi
  for i=1:4
    for j=1:4
      for c=1:4
        tmp(:,:,:,i,j) = tmp(:,:,:,i,j) + gt(:,:,:,i,c).*gi(:,:,:,c,j);
      end
    end
  end

  for i=1:4
    for j=1:4
      for c=1:4
        git(:,:,:,i,j) = git(:,:,:,i,j) - gi(:,:,:,i,c).*tmp(:,:,:,c,j);
      end
    end
  end

end
